Supplementary Fig 8 - GEP usage per thymic lineage
Salomé Carcy
30 July 2024
1 IMPORT
1.1 Import librairies
library(ggplot2)
library(cowplot)
library(tidyverse)
library(dplyr)
library(Seurat)
setwd("~/Projects/HumanThymusProject/")
source("~/Projects/HumanThymusProject/scripts/colors_universal.R")1.2 Import data
# import seurat objects
seur.nkt <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_inkt.rds")
seur.mait <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_mait.rds")
seur.gdt <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_gdt.rds")
seur.cd4 <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd4.rds")
seur.cd8 <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd8.rds")
# sanity check seurat objects correspond to Fig 2
DimPlot(seur.nkt, group.by="clusters_per_lineage", label=T, cols=cols_thym_nkt)
DimPlot(seur.mait, group.by="clusters_per_lineage", label=T, cols=cols_thym_mait)
DimPlot(seur.gdt, group.by="clusters_per_lineage", label=T, cols=cols_thym_gdt)
# import GEP usage dataframe
gep_usage <- read.table("./data_github/cNMF/non_imputed_cNMF_allcells.usages.k_12.dt_0_02.consensus.txt", header=T)
colnames(gep_usage) <- paste0("GEP", c(2,5,3,1,4,12,6,7,8,10,9,11), "_usage")
gep_usage <- gep_usage[,paste0("GEP", 1:12, "_usage")]
head(gep_usage)## GEP1_usage GEP2_usage GEP3_usage GEP4_usage GEP5_usage
## CD4_1_Thymus_878524 0.0000000 0.00000000 0.01697069 0.02701780 0.003293012
## CD4_1_Thymus_679405 0.0000000 0.00000000 0.00000000 0.00000000 0.000000000
## CD4_1_Thymus_281188 0.0000000 0.01499587 0.02538062 0.01897941 0.000000000
## CD4_1_Thymus_285601 0.3893558 0.19701648 0.01374016 0.03301096 0.009603295
## CD4_1_Thymus_174537 0.0000000 0.00000000 0.01813514 0.02750585 0.001959355
## CD4_1_Thymus_128406 0.0000000 0.02029776 0.02094953 0.01324090 0.001427132
## GEP6_usage GEP7_usage GEP8_usage GEP9_usage GEP10_usage
## CD4_1_Thymus_878524 0 0.0000000000 0.01670451 0.06574529 0.14537527
## CD4_1_Thymus_679405 0 0.0000000000 0.00000000 0.00000000 0.00000000
## CD4_1_Thymus_281188 0 0.0002823688 0.01075830 0.09749147 0.21138272
## CD4_1_Thymus_285601 0 0.0178349300 0.22377394 0.05961335 0.01831705
## CD4_1_Thymus_174537 0 0.0000000000 0.01610514 0.06646363 0.15172973
## CD4_1_Thymus_128406 0 0.0055636405 0.01326269 0.08876818 0.21684438
## GEP11_usage GEP12_usage
## CD4_1_Thymus_878524 0.655783300 0.00000000
## CD4_1_Thymus_679405 1.147190900 0.00000000
## CD4_1_Thymus_281188 0.563456800 0.00000000
## CD4_1_Thymus_285601 0.008145613 0.03549316
## CD4_1_Thymus_174537 0.649509700 0.00000000
## CD4_1_Thymus_128406 0.565640400 0.000000002 ANALYSIS
2.1 Add GEP1-GEP11 usage to seurat objects
# check that rownames are in same order
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data),])==rownames(seur.cd4@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data),])==rownames(seur.cd8@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data),])==rownames(seur.nkt@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data),])==rownames(seur.mait@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data),])==rownames(seur.gdt@meta.data), useNA="ifany") # all TRUE
# add GEP usage to seurat objects
gep_colnames <- paste0("GEP", 1:11, "_usage")
seur.cd4@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data), gep_colnames]
seur.cd8@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data), gep_colnames]
seur.nkt@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data), gep_colnames]
seur.mait@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data), gep_colnames]
seur.gdt@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data), gep_colnames]2.2 Plot GEP usage per thymic lineage
# Plot
seur_vector <- list("iNKT"=seur.nkt, "MAIT"=seur.mait, "GDT"=seur.gdt, "CD4"=seur.cd4, "CD8"=seur.cd8)
plist_massive <- list()
for(gep in gep_colnames){
# print(gep)
# get row of plots (for one GEP)
plist <- list()
for(i_seur in names(seur_vector)){
# print(i_seur)
seur <- seur_vector[[i_seur]]
# get legend once
if(gep=="GEP1_usage" & i_seur=="iNKT"){
plegend <- ggpubr::get_legend(
SCpubr::do_FeaturePlot(seur, features=gep, ncol=1, legend.position="right", legend.title="GEP usage", border.size=1, pt.size=2, order=T)+
scale_color_viridis_c(limits=c(0,1.25), option="B")+
theme(plot.background = element_rect(color = "black"))
)
}
# make plot
legendpos <- "none"
plt_title <- ""
if(gep=="GEP1_usage"){plt_title=i_seur}
p <- SCpubr::do_FeaturePlot(seur, features=gep, ncol=1, legend.position=legendpos, border.size=1, pt.size=2, order=T, plot.title=plt_title)+
scale_color_viridis_c(limits=c(0,1.25), option="B")+
theme(plot.background = element_rect(color = "black", linewidth = 2),
plot.margin=unit(c(0,5.5,0,5.5), "points"))
plist[[i_seur]] <- ggrastr::rasterise(p, layers="Point", dpi=300)
}
prow_gep <- plot_grid(plotlist=plist, nrow=1, scale=0.9)
plist_massive[[gep]] <- prow_gep
}
plot_grid(
plot_grid(plotlist=plist_massive, ncol=1, align="h"),
ggpubr::as_ggplot(plegend), ncol=2, rel_widths = c(5, .5), scale=c(0.95, 5)
)
3 SESSION INFO
sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] SeuratObject_4.1.3 Seurat_4.3.0.1 forcats_1.0.0 stringr_1.5.0
## [5] dplyr_1.1.2 purrr_1.0.1 readr_2.1.4 tidyr_1.3.0
## [9] tibble_3.2.1 tidyverse_1.3.2 cowplot_1.1.1 ggplot2_3.4.2
##
## loaded via a namespace (and not attached):
## [1] readxl_1.4.2 backports_1.4.1 plyr_1.8.8
## [4] igraph_1.5.0 lazyeval_0.2.2 sp_2.0-0
## [7] splines_4.1.3 listenv_0.9.0 scattermore_1.2
## [10] digest_0.6.32 htmltools_0.5.5 viridis_0.6.3
## [13] fansi_1.0.4 magrittr_2.0.3 tensor_1.5
## [16] googlesheets4_1.1.1 cluster_2.1.4 ROCR_1.0-11
## [19] tzdb_0.4.0 globals_0.16.2 modelr_0.1.11
## [22] matrixStats_1.0.0 timechange_0.2.0 spatstat.sparse_3.0-2
## [25] colorspace_2.1-0 rvest_1.0.3 ggrepel_0.9.3
## [28] haven_2.5.3 xfun_0.39 crayon_1.5.2
## [31] jsonlite_1.8.7 progressr_0.13.0 spatstat.data_3.0-1
## [34] survival_3.5-5 zoo_1.8-12 glue_1.6.2
## [37] polyclip_1.10-4 gtable_0.3.3 gargle_1.5.1
## [40] leiden_0.4.3 car_3.1-2 future.apply_1.11.0
## [43] abind_1.4-5 scales_1.2.1 DBI_1.1.3
## [46] SCpubr_1.1.2 rstatix_0.7.2 spatstat.random_3.1-5
## [49] miniUI_0.1.1.1 Rcpp_1.0.10 viridisLite_0.4.2
## [52] xtable_1.8-4 reticulate_1.30 htmlwidgets_1.6.2
## [55] httr_1.4.6 RColorBrewer_1.1-3 ellipsis_0.3.2
## [58] ica_1.0-3 pkgconfig_2.0.3 farver_2.1.1
## [61] sass_0.4.6 uwot_0.1.16 dbplyr_2.3.2
## [64] deldir_1.0-9 utf8_1.2.3 tidyselect_1.2.0
## [67] labeling_0.4.2 rlang_1.1.1 reshape2_1.4.4
## [70] later_1.3.1 munsell_0.5.0 cellranger_1.1.0
## [73] tools_4.1.3 cachem_1.0.8 cli_3.6.1
## [76] generics_0.1.3 broom_1.0.5 ggridges_0.5.4
## [79] evaluate_0.21 fastmap_1.1.1 yaml_2.3.7
## [82] goftest_1.2-3 knitr_1.43 fs_1.6.2
## [85] fitdistrplus_1.1-11 RANN_2.6.1 pbapply_1.7-2
## [88] future_1.33.0 nlme_3.1-162 mime_0.12
## [91] ggrastr_1.0.2 xml2_1.3.4 compiler_4.1.3
## [94] rstudioapi_0.14 beeswarm_0.4.0 plotly_4.10.2
## [97] png_0.1-8 ggsignif_0.6.4 spatstat.utils_3.0-3
## [100] reprex_2.0.2 bslib_0.5.0 stringi_1.7.12
## [103] highr_0.10 lattice_0.21-8 Matrix_1.5-4.1
## [106] vctrs_0.6.3 pillar_1.9.0 lifecycle_1.0.3
## [109] spatstat.geom_3.2-1 lmtest_0.9-40 jquerylib_0.1.4
## [112] RcppAnnoy_0.0.21 data.table_1.14.8 irlba_2.3.5.1
## [115] httpuv_1.6.11 patchwork_1.1.2 R6_2.5.1
## [118] promises_1.2.0.1 KernSmooth_2.23-21 gridExtra_2.3
## [121] vipor_0.4.5 parallelly_1.36.0 codetools_0.2-19
## [124] assertthat_0.2.1 MASS_7.3-60 withr_2.5.0
## [127] sctransform_0.3.5 parallel_4.1.3 hms_1.1.3
## [130] grid_4.1.3 rmarkdown_2.23 carData_3.0-5
## [133] googledrive_2.1.1 Cairo_1.6-0 Rtsne_0.16
## [136] ggpubr_0.6.0 spatstat.explore_3.2-1 shiny_1.7.4
## [139] lubridate_1.9.2 ggbeeswarm_0.7.2